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In the framework of various statistical models as well as of mechanisms for color confinement, 
disorder parameters can be developed which are generally expressed as ratios of partition functions 
and whose numerical determination is usually challenging. We develop an efficient method for their 
computation and apply it to the study of dual superconductivity in Ad compact f/(l) gauge theory. 
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I. INTRODUCTION 

Order-disorder transitions are common to a wide class 
of models in Statistical Mechanics and Quantum Field 
Theory, the Ising model being a prototype tli,^2]. In those 
models, one phase is characterized by the condensation of 
dual topological excitations which spontaneously breaks 
a dual symmetry, and correlation functions of those ex- 
citations can serve as disorder parameters for the transi- 
tion: in general they are non local in terms of the original 
variables, so that their numerical study can be challeng- 
ing. The difficulty becomes evident when the correlation 
functions are expressed as ratios of partition functions. 

Relevant examples are encountered when studying 
color confinement in QCD: that is usually believed to be 
related to the condensation of some topological excita- 
tions and models can be constructed accordingly, which 
place the confinement-deconfinement transition into the 
more general scenario of order-disorder transitions. One 
appealing model is based on dual superconductivity of 
the QCD vacuum and relates confinement to the breaking 
of an abelian dual symmetry induced by the condensation 
of magnetic monopoles [1, 0, IB ■ The possibility to de- 
fine disorder parameters in this scenario has been studied 
since a long time [1, 0, H, • One parameter has been 
developed by the Pisa group and is the expectation value 
of an operator which creates a magnetic monopole; 
(^) has been shown to be a good parameter for confine- 
ment in U{1) [Toll in pu re Yang-Mills theories [ll|, ^ 
and in full QCD |l3l . [l4| : similar parameters have been 
developed both in gauge theories [H, [H, EBl and in 
statistical models [H, US IHI ■ The operator ^ is express- 
ible as the exponential of the integral over a time slice, 
/i = exp(— /^AS") (see later for details), so that its expec- 
tation value can be rewritten (we consider a pure gauge 
theory as an example) as 



,-/3(S+AS) 



(1) 



where the functional integration is over the gauge link 
variables, S is the euclidean action of the theory and /3 
is the inverse gauge coupling. The difficulty involved in 



its numerical computation stems from the poor overlap 
among the two statistical distributions corresponding to 
the partition functions Z and Z: configurations which 
give significant contributions to Z are instead extremely 
rare in the original ensemble corresponding to Z, so that 
they are very badly sampled in a Monte Carlo simulation. 
The problem worsens rapidly when increasing the spatial 
volume, making a determination of (fi) hardly feasible. 
One way out is to evaluate susceptibilities of (/i), like: 



d(3 



ln(M) - {S)s - {S)s , 



(2) 



from which the disorder parameter can eventually be re- 
constructed as follows 



(m> exp 



p(/3')d/?' 



(3) 



While that is enough to test (/i) as a parameter for con- 
finement, a direct determination could be useful in con- 
texts like the study of its correlation functions [l^, [H, 

US- 

The problem of dealing with extremely rare configu- 
rations can be app roached using the idea of generalized 
ensembles [2^ 113 Izl]- In that framework we propose 
a new method for a direct computation of (fx), which is 
inspired by analogous techniques used for the study of 
the 't Hooft loop [22|. We describe the method for the 
case of the 4d compact U{1) gauge theory in the Wilson 
formulation, but it is applicable to the study of disorder 
parameters in a wide class of analogous problems. 



II. THE METHOD 

The partition function of the model is defined, in the 
Wilson formulation, as follows 



Z(/3)- /He-^^ 



S= (l-cos(V(5',i)) 



(4) 



(5) 
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where the integration is over the Hnk variables (phases in 
U{1)) and 9^^{x,t) is the plaquette in the fiv plane sit- 
ting at lattice site {x,t). The model has a critical point 
at — 1.01, which is believed to be weak first order and 
separates a disordered phase (/3 < [3c) ^ with condensa- 
tion of magnetic monopoles and confinement of electric 
charges, from a Coulomb phase where magnetic charge 
condensation disappears. 

The magnetically charged operator ^(y^to), whose ex- 
pectation value detects dual superconductivity, creates a 
monopole in y at time to by shifting the quantum gauge 
fields by the classical vector potential of a monopole, 
b{x — y), and can be written (see Ref. [l3] for details) 
as 



ti-{y,to) = exp 



.1 



d^x E{x, to)h{x — y) 



(6) 



with the electric field E{x, to) being the momentum con- 
jugate to the quantum vector potential. It can be dis- 
cretized on the lattice as follows: 

^ ^ g/3Ex,,(cos(eoi(5,to)-6.(5-tr))-cos(9oi(£,to))) ^ e"''^'^ (7) 

where 9oi are the phases of the temporal plaquettes, cor- 
responding to the electric field in the (nai've) continuum 
limit. If we define the modified action 

S= (1 - cos(0^.)) 

itMto,(Mi^) 

+ ^ (1 - cos(%(f , to) ~ b,{x - y))) ^ S + AS (8) 



which differs from S only on the timeslice where the 
monopole has been created, we can write 



{^Ji{y,t)) = Z{^3)/Z{P)■ 
Z{p) = / [dO] exp(-/3S') . 



(9) 
(10) 



Measuring (/x(y, t)) in a Monte Carlo (MC) simula- 
tion is very difficult [lO|: fi gets significant contributions 
only on those configurations having very small statistical 
weight (which are poorly sampled in a finite MC simu- 
lation). The difhculty increases with the system size as 
the two distributions corresponding to Z and Z shrink 
towards non overlapping delta functions in the configu- 
ration space. Our proposal is to determine the ratio in 
Eq. ([9]) by using intermediate distribution functions hav- 
ing a reasonable overlap with both statistical ensembles 
corresponding to Z and Z: our method has many sim- 
ilarities with strategies adopted in the computation of 
analogous order parameters [1^, [2^, [13, HH . As a first 
step we rewrite the ratio as the product of distinct 
ratios: 



Z _ Zn Zn-1 Z\ 
Z Zn-1 Zn^2 Zo 



where Zm = Z ^ Zo = Z and Z^ is defined in terms of an 
action Sj. which is an interpolation between S and S: 



Zk = J[d0]t 



9 -^-^-^ ^ 



(12) 
(13) 



The idea is to compute each single ratio by a different 
Monte Carlo simulation: the difhculty of dealing with 
N simulations should be greatly compensated by the in- 
creased overlap in the distributions corresponding to each 
couple of partition functions, leading to a benefit which 
increases exponentially with N. As a second step to fur- 
ther improve the overlap, we compute each single ratio on 
the r.h.s. of Eq. (jlip using an intermediate distribution: 



■'k+l 



(exp(-/3A5/27V)),+i/2 



(14) 



Zk (exp(/3A5/2A^)),+i/2 ' 
where each expectation value is computed with the action 
Sk+i/2 = (1 - (fc + l/2)/A^) S+{{k + 1/2) IN) S. (15) 



(11) 



Since both expectation values in Eq. (jl4p are computed 
with the same MC simulation, we make use of a jack- 
knife analysis to get a reliable error on Zk+i/Zk- The 
final uncertainty on Z /Z is then obtained by standard 
error propagation since each single ratio on the r.h.s. of 
Eq. (fTT|) is obtained by an independent MC simulation. 

Our technique of rewriting the ratio Z /Z as a prod- 
uct of intermediate ratios resembles very closely the well 
known snake algorithm [soj as well as other algorithms 
inspired by it, like that used for the computation of the 
helicity modulus [3l| . However it differs from previous al- 
gorithms in the choice of the intermediate partition func- 
tions, which in our case is not related to the details of the 
model, so that it can be applied without modifications 
to a wider class of problems in Lattice Field Theory and 
Statistical Mechanics. 

Regarding the choice of boundary conditions (b.c), we 
do it in a consistent equal way for all the partition func- 
tions in Eq. pT|) . We make use of both periodic and free 
b.c. in the spatial directions: while one could expect a 
substantial difference in presence of a magnetic charge, 
we will show that, in the phase where (/i) ^ 0, the two 
choices lead to the same thermodynamical limit; that is 
expected since in that phase the vacuum does not have 
a well defined magnetic charge and a monopole is com- 
pletely screened. As for the temporal direction, a consis- 
tent usual choice for (/i) [l3| is that of periodic b.c. for Z 
and C* b.c. for Z; in particular C* boundary conditions, 
which corresponds to performing a charge conjugation 
transformation on gauge fields when crossing the time 
boundary, are taken so as to annihilate the monopole 
after one loop around the periodic time direction, avoid- 
ing in this way that it propagates an indefinite number of 
times. However we do not keep that choice, since it would 
lead to intermediate actions with inconsistent mixed b.c; 
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instead we adopt either free or periodic b.c. for both Z 
and Z also in the time direction, showing again that the 
choice is inessential when (/i) 7^ 0. 

III. NUMERICAL RESULTS 

As a first test we compare the nai've computation of 
(/i), i.e. performed with a single MC simulation using the 
Wilson action S*, with our method for = 1: we use a 
4'* lattice with free b.c. at /3 = 0.8 and 10^ measurements 
for both cases, obtaining (/i) = 1.14(18) with the naive 
computation and (/i) — 0.868(3) with our method. Apart 
from the strongly reduced error, much is learned by look- 
ing at the distributions of the observables exp(— /3AS')s 
and exp(±/3AS'/2)^g_^j^^2 (the subscript indicates the ac- 
tion used for sampling) used in the computation (see 
Eg. fH)) . In Fig.[l]we plot, for each fixed observable O, the 
distribution of the logarithm of O times the observable 
itself (we choose the logarithm for graphical convenience) 
as a function of log O, so that the integral of each curve 
gives the expectation value of the relative observable. As 
it is clear, for (exp(— /3AS'))5 most of the contribution 
comes from a region which is badly sampled: on larger 
lattices the problem worsens rapidly and a naive deter- 
mination of (/i) is unfeasible. The improvement obtained 
with our method is apparent already for A^ = 1. 

In Fig. [2] we show a determination of (/i) for several 
values of A^ on a lattice 16'' with free b.c. at /3 = 0.8: 
for each determination a comparable whole statistics of 
A^ X A'mcas = 3.2 • 10^ measurements has been used, so 
that the error on (/i) is an indication of the efficiency 
as a function of A^; in Fig. [3] we report the intermediate 
ratios Zk+i/Zk (see Eq. [11]) used for each measurement. 
While the intermediate ratios are strongly dependent on 
A^, (/J,) is not, thus confirming the absence of uncontrolled 
systematic errors. The statistical error rapidly changes 
for small values of A^, but then stabilizes, indicating that 
a value A^ ~ 0(10) saturates the improvement. 

As an application of our method we analyse some rel- 
evant features of the disorder parameter, starting with a 
study of its thcrmodynamical limit in the confined phase. 
In Fig. d] we show {ji) determined with both free and pe- 
riodic b.c. at /3 = 0.5 as a function of the lattice size L. 
A fit according to {^) ^ A + B/L gives A = 0.945(6) and 
B = 0.31(6) with free b.c. (x^ = 0.5) and A = 0.940(6) 
and B = 0.2(6) with periodic b.c. (x^ = 0.8). In both 
cases (n) has a well defined thermodynamical limit, which 
does not depend (within numerical errors) on the b.c. 
chosen: that is expected in the phase where magnetic 
charge is completely screened. We stress that, contrary 
to what may happen with other parameters 31], we do 
not expect exactly (^) = 1 in the confined phase: indeed 
any non-zero value of the disorder parameter ensures the 
breaking of the magnetic symmetry, hence dual super- 
conductivity. 

A further confirmation of screening comes from the 
study of cluster property in the correlation functions. In 
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FIG. 1: For each observable O involved in the computation 
on the 4* lattice at /3 = 0.8, we plot O times the distribution 
of log O as a function of log O: the integral under each curve 
gives (O). 

Table |T] we report the values measured for various tem- 
poral and spatial correlators of /i using periodic b.c, to- 
gether with their second and fourth powers. In this way 
we compare, for instance, the value of (/i)^ (first row, sec- 
ond column) with that of the two point function (/I/i) at 
large distances (first column, second row), or (yit)'* with 
the four point function, and so on: the compared quan- 
tities should approach each other (exponentially in the 
extension of the higher order correlator) if cluster prop- 
erty is obeyed. This is nicely verified, within errors, from 
the data reported in the Table. 

Results are quite different in the deconfined phase. In 
Fig. [5] we show {^) as a function of i at /3 = 1.1: the 
determinations with free and periodic b.c. differ from 
each other, both going to zero exponentially with the 
lattice size L. This is the correct expected behavior in 
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FIG. 2: Determination of {^) on a lattice 16* with free b.c. 
at /3 = 0.8 as a function of N . 
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FIG. 3: Intermediate ratios Z^+i/Zk used for the determina- 
tions reported in Fig. [5] 
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0.439(12) 0.193(11) 0.037(4) 

0.182(7) 0.033(3) 

0.183(12) 0.033(4) 
0.037(6) 



TABLE I: Determination a,t (3 — 0.8 on a 16* lattice of (fj), 
of its spatial and temporal 2-point function (second and third 
row) and of its mixed 4-point function (last row) , with t = 8 
and z= (0,0,8). The measured correlator is reported in the 
first column, indicated generically with O, while in the sec- 
ond and third column (indicated with and C*) we report 
respectively the second and fourth power of the correlator, 
which are used to test cluster property on the higher order 
correlators reported in lower rows. Within errors all measure- 
ments are compatible with the hypothesis that the correlators 
are already in their asymptotic regime governed by cluster 
property. 



IV. CONCLUSIONS 



the phase with magnetic charge superselection [l^, [S^l . 

Finally we consider the critical behaviour of the dis- 
order parameter close to the phase transition, where 
(/i) ~ (r is the reduced temperature). That trans- 
lates in the following finite size scaling (f.s.s.) behaviour 



(16) 



where </) is a scaling function. To test this ansatz we have 
determined (/i) close to the phase transition on several 
different lattice sizes. Fixing the known value of (3c = 
1.011 and u = 1/d ~ 0.25 as appropriate for a weak 
first order transition, we obtain a reasonable scaling with 
S 2.3: the quality of our f.s.s. analysis is shown in 
Fig.O 



We have proposed a new technique for the computation 
of disorder parameters and applied it to the study of the 
parameter for dual superconductivity in 4ci compact U(l) 
gauge theory. Our method is inspired by methods used 
for the study of the 't Hooft loop [soj . 

We have determined some relevant features of (fi) both 
in the confined and in the Coulomb phase. A careful 
analysis of its critical properties could help in clarify- 
ing the nature of the pha se transition at zero as well as 
at finite temperature j33l . [s^ : to that aim also a direct 
comparison with analogous order parameters developed 
for C/(l) jl8, 31] will be particularly useful. 

Our method can be placed in the more general frame- 
work of techniques based on the idea of generalized en- 
sembles 23, 24, 25]: in that respect, it has the advantage 
to provide a recipe which can be easily applied, with 
none or few modifications, to a wide class of problems. 
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FIG. 4: (^) a,t /3 — 0.5 as a function of the lattice size L 
determined with both free and periodic b.c. 



FIG. 5: {n) at /3 = 1.1 as a function of the lattice size L 
determined with both free and periodic b.c. 
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FIG. 6: F.s.s. analysis of (fi) around the phase transition 



Among others we will consider in the future the study of 
order-disorder transitions in statistical models and dual 



superconductivity in non Abelian gauge theories. An- 
other benefit of our proposal is that it leaves room for 
considerable further improvement: for instance it could 
be possible to choose a more general non linear interpo- 
lation between the two actions S and S, differently from 
what has been done in Eq. ([T3|) , so as to concentrate the 
numerical effort on those intermediate ensembles where 
the statistical distribution is changing more rapidly. 
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